NASA 

Technical 

Paper 

2919 


1989 


Interactions of 
Tollmien-Schlichting 
Waves and Dean 
Vortices 

Comparison of Direct Numerical 
Simulation and a Weakly 
Nonlinear Theory 


Bart A. Singer 

High Technology Corporation 

Hampton , Virginia 

Thomas A. Zang 
Langley Research Center 
Hampton , Virginia 


NASA 

National Aeronautics and 
Space Administration 

Office of Management 

Scientific and Technical 
Information Division 


1. Introduction 


The study of incompressible curved channel flow can illuminate some important issues dealing 
with laminar-turbulent transition. The streamwise curvature of the flow induces an instability 
resulting in longitudinal vortices (hereafter referred to as Dean vortices (ref. 1 )). For mildly curved 
channels, the minimum Reynolds number for instability of the Dean vortices is comparable with 
that for Tollmien-Schlichting (TS) waves. In this regime, Gibson and Cook (ref. 2 ) showed that 
oblique waves are never the dominant linear disturbance. This makes the flow an ideal candidate 
for studying the nonlinear interactions between two-dimensional ( 2 D) TS waves and the streamwise- 
oriented Dean vortices. The fact that the flow is strictly parallel implies that a self-consistent analysis 
can be performed without the application of triple-deck theory (ref. 3). 

To date, the only weakly nonlinear theory for the interaction of Dean vortices and TS waves 
at finite Reynolds numbers is that developed by Daudpota, Hall, and Zang (hereafter referred to 
as DHZ) (ref. 3). They employ a multiple-scale version of the Stuart (ref. 4) and Watson (ref. 5) 
approach to approximate the complex equations describing the flow as two coupled Landau equations 
for the perturbation amplitudes of Dean vortices and TS waves. These Landau equations are 
examined in the context of a nonlinear autonomous system (for details on nonlinear autonomous 
systems see, for example, ref. 6 ) to determine the characteristics of the various equilibrium states 
(nodes). Using this information, one can make predictions about the behavior of the real flow for 
ranges of parameters. The limits of these ranges are not known a priori; they must be determined 
by some other means. Direct numerical simulation (DNS) of the full equations describing the flow 
can be used to check the validity of the theory at specific points in the parameter space. In this 
paper we use DNS to determine the accuracy of the weakly nonlinear theory developed by DHZ. 

In section 2 we enumerate the relevant nomenclature and describe the weakly nonlinear theory. 
Section 3 describes the DNS code with which we compare the weakly nonlinear predictions. We 
compare the results in section 4 and discuss their implications in section 5. 

The authors would like to thank Q. Isa Daudpota and Philip Hall for their many useful discussions 
associated with this project. 

2. Mathematical Formulation of the Weakly Nonlinear Theory 

Here we introduce the notation and scalings used in the weakly nonlinear theory. The notation 
is the same as that used by DHZ. The interested reader is referred there for a detailed description 
of the theory; here we give only the information that is essential for its application. 

Cylindrical coordinates (r, 6 , z) are used to describe the flow between infinite concentric walls of 
outer radius R 0 and inner radius R Figure 1 illustrates the geometry of the problem. The distance 
between the walls d = R 0 — Ri is used to nondimensionalize the axial (spanwise) coordinate 2 , while 
the outer radius R 0 nondimensionalizes the radial coordinate r. The ratio 77 = Ri/R 0 is used to 
describe the curvature of the channel; the radial coordinate r varies between 77 and 1. The quantity 
1 - 77 appears often and is denoted by A. At the critical curvature parameter X c = 2.179 x 10“ 5 , 
the minimum Reynolds number for instability of TS waves is the same as that for Dean vortices 
(see Gibson and Cook, ref. 2). When plotting velocity distributions in the radial direction, it is 
convenient for the abscissa to vary between —1 and 1 , hence we define another radial coordinate, 
y = (2r - (1 + y)) /A. 

The mean flow is denoted by radial velocity U, azimuthal velocity V, axial velocity W, and 
pressure P. When a constant azimuthal pressure gradient ^ drives the flow, a solution to the 
momentum equations gives (U,V,W) as 
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The kinematic viscosity v and the density p are constant. The Reynolds number is defined as 
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(4) 


In the limit of a straight channel, V m becomes four times the mean centerline velocity, and the 
Reynolds number approaches eight times the standard straight channel Reynolds number based on 
centerline velocity and half-width. 

The mean flow described above is disturbed with pressure p and velocity (u, v, w) perturbations. 
The perturbation pressure is scaled with pV^, the radial and axial velocities with the diffusive scale 
v/d, and the azimuthal velocity v with the convective velocity scale V m . Time t is normalized with 
Ro/'Vm- 

On the neutral stability curve the linear perturbations due to the Dean vortex are proportional 
to eexp(ikz), where e is an arbitrary small parameter, i = 1 , and the spanwise wave number k is 

scaled with channel width d . Linear TS disturbances are proportional to e exp(im0) exp(crt), where 
the nondimensional complex growth rate a is defined as 


a — °R + 


(5) 


and (tr = 0 on the TS neutral stability curve. The dimensional streamwise wave number for a TS 
wave, ajSi re lated to the wave number m via the relation 


a^sdx = mdO 


( 6 ) 


where dx is the differential distance along the centerline of the channel. Hence, using 

dx = — + — dO (7) 

one obtains the relation 

d ' Q TS = (8) 

where the quantity d is included so that Qxs has the appropriate dimensions. Note that d • a^s is 
nondimensional such that its value is twice the standard straight channel value of the nondimensional 
streamwise wave number. 

Neutral stability curves for TS and Dean disturbances for curvatures slightly less than, equal 
to, and slightly greater than the critical value are shown in figures 2(a) and (b). The apparent 
sensitivity of the TS disturbances to curvature is an artifact of using m as the abscissa. Using the 
transformation equation (8) one finds that the plots overlap; differences in the value of the neutral 
Reynolds number for a given TS wavelength are less than 20 for the three cases shown. For the 
Dean disturbances, differences in the neutral Reynolds number are about 200. For consistency with 
standard nomenclature, the left branch of these curves will be called the lower branch and the right 
side will be called the upper branch. Reference 3 concentrates on the lower branch of the TS neutral 

1 Typographical errors in DHZ account for a missing factor of 77 in their definition of f(r). 
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curve because it is the more relevant branch when one considers boundary-layer-type flows. In the 
case of the Dean instability, harmonics of the lower branch wave number can fall in the unstable 
region and actually grow faster than the primary disturbance considered by the weakly nonlinear 
theory. To avoid this possibility, here we restrict our attention to the Dean upper branch. 

The Reynolds number is expanded about a value Rq that is on both the TS and the Dean neutral 
stability curves so that 

Re = Rq + R\ 

The curvature parameter A is written as a deviation from the critical curvature ratio at which the 
minimum Reynolds number is the same for TS and Dean disturbances. Hence 


A = A c + e 2 A 


A slow time variable is defined by 

r = eh 

The quantity eX(r) represents the time history of half the maximum azimuthal velocity of the Dean 
vortex, evaluated on the channel centerline. The corresponding quantity for TS waves is eY(r), 
which measures half the maximum radial velocity on the centerline. The following coupled set of 
Landau equations were derived by DHZ: 

\ = (Rift + All) m 2 + h\X\ 2 \X \ 2 + Vl\X\ 2 \Y\ 2 (9) 

\ = (Rilh + A72) l*f + Wl 2 !*! 2 + m\Y\ 2 \v \ 2 (10) 

The ratios — 7i//?i and — 72/^2 are the slopes of the zero-growth-rate curves for \X\ and |F| in 

the (A,i?i) plane. Numerical computations (confirmed by us) indicate that A 72 << R 1 P 2 for 

all cases considered, so 72 is taken to be zero in what follows. The seven remaining coefficients, 
71 , A, /? 2 i #i, 62 , r/i, and 77 2 , are all computed at points where the Dean and TS neutral stability 
curves intersect. These values were computed by DHZ, using 51 points across the channel with 
a fourth-order finite-difference discretization. We used the same software and recomputed the 
coefficients, using 501 points across the channel. All but the coefficient 7/2 remained constant to 
three digits with the increased resolution; 772 changed by about 20%. 

Equations (9) and (10) admit four possible steady state solutions: the trivial solution, finite \X\ 
with \Y\ = 0, finite \Y\ with \X\ = 0, and a combined state with finite values of both \X\ and \Y\. 
Each of these states is a critical point and may be classified according to the local behavior of the 
solutions in the vicinity of the point. In the phase plane (\X\ - |y|), stable nodes are identified as 
those for which solution trajectories near the point converge to the node, unstable nodes are critical 
points from which solution trajectories diverge, and saddle points are those points for which a finite 
number of trajectories converge to the point while all others diverge. For the case of a channel more 
strongly curved than critical (A > 0), DHZ found that the trivial solution is an unstable node, the 
TS equilibrium state is a saddle point, and the Dean equilibrium state is a stable node. No combined 
mode solution exists here. For any initial condition that contains even a very small amplitude Dean 
disturbance, the solution will ultimately tend towards an equilibrium Dean disturbance with no TS 
perturbation. The situation with A < 0 is more complicated. For small Reynolds number deviations 
from the neutral curve, only the TS equilibrium state will exist. At larger deviations, the TS state 
itself is unstable and the combined state of TS and Dean disturbances is stable. For sufficiently large 
Reynolds number perturbations, only the equilibrium Dean disturbance will persist. 

3. The Direct Numerical Simulation 

Here we describe the numerical techniques used to solve the full Navier-Stokes equations in 
the curved channel. The basic equations are written in perturbation variables; the total velocity is 
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(U, V \ W) + ( u *, v*, iu*), where ([/, V, VF) are given by equations (1)— (3) of section 2. The coordinates 
and velocities are then nondimensionalized with the outer radius of curvature R 0 and the quantity 
V m , which is defined in equation (2). Primed variables are the nondimensional counterparts of the 
respective starred quantities. In primitive variable form, the continuity and momentum equations 
are 
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where 

v 2 — { —} — — — 

r dr\dr) + r 2 86 2 dz 2 


(15) 


The quantity / is the same as that defined in equation (3) and represents the nondimensional mean 
flow. It is computed using double precision to ensure sufficient accuracy. The Reynolds number 
used here is the same as that defined by equation (4). The relationship between (u', v\ w f ,p f ) and 
(u, w,p) of section 2 is 


u = uf Re 
v ! = v 
w — w/Re 

P = P (16) 

The boundary conditions are periodic in the azimuthal and axial directions. All velocity components 
are required to be zero at the inner and outer channel walls. 

To solve these equations, we used a curved channel variant of the program described by Zang 
and Hussaini (refs. 7 and 8). It employs Fourier expansions in the 9 and z directions and Chebyshev 
expansions in r. The independent variables are advanced in time by using Crank-Nicholson for the 
diffusion, pressure gradient, and incompressibility constraint terms, and third-order Runge-Kutta 
for the remaining terms. The CFL number limitation due to the stability of the Runge-Kutta 
method is y/3/ir. The pressure is computed on a staggered grid, and an advanced splitting method 
described by Zang and Hussaini (ref. 8) is used to reduce the errors at the walls. The methods 
require the solution of a Poisson equation for the pressure, and Helmholtz equations for the velocity 
components. Typically, at the end of a full time step, either the incompressibility constraint or 
the no-slip boundary condition is satisfied, but not both. However, our code enforces boundary 
conditions on the intermediate steps such that the error due to splitting the equations is reduced to 
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O ^A£ 2 J. (See Zang and Hussaini (ref. 8) for details.) The Helmholtz equations are solved with a 
Chebyshev tau method and the pressure equation by an iterative approach. 

The spatial derivatives are approximated with analytic derivatives of the series representations. 
The nonlinear terms are evaluated most efficiently in physical space. Fast Fourier and Chebyshev 
transforms permit the data to be rapidly converted to and from its physical and spectral space 
representations. In this approach, as with other collocation methods, both truncation and aliasing 
errors are present. Canuto, Hussaini, Quarteroni, and Zang (ref. 9) discuss both types of errors 
in detail and provide a thorough discussion of available spectral algorithms for channel flow 
simulations. Recently, Zang (ref. 10) demonstrated that numerical simulations using collocation 
schemes experience smaller aliasing errors when the nonlinear convective terms are written in the 
skew-symmetric form rather than the rotation form. In the simulations discussed here, the skew- 
symmetric form has been used. On a 64 3 grid, the code uses approximately 90 /isec per grid point 
per step on a single Cray-2 processor. 

The DNS code was extensively tested before being used for this project. Linear eigenmodes for 
both Dean and TS disturbances were obtained from a linear stability code used by DHZ. This code 
employs 501 grid points in a finite-difference method developed by Malik, Chuang, and Hussaini 
(ref. 11). At the critical curvature parameter A c = 2.179 x 10“ 5 , 2D TS complex growth rates 
from the curved channel linear stability code differ from the corresponding quantities obtained 
from a straight channel spectral Orr-Sommerfeld solver by approximately 0.1%. A comparison of 
instantaneous complex growth rates from the DNS code with those obtained by the curved channel 
stability code is shown in table I. All cases were performed at the critical curvature ratio. The DNS 
data were obtained by using 48* Chebyshev modes in the radial direction. Initial amplitudes were 
1.0 x 10 -8 of the mean flow centerline velocity. The growth rates and frequencies were obtained 
after 100 time steps with At = 10~ 5 . Differences between the linear theory and the DNS are more 
pronounced for the TS waves (m ^ 0) than for the Dean disturbances (m = 0), but in all cases the 
differences are within the uncertainty range of the linear code. 

Table I. Complex Growth Rates 


Re 

m 

k 

Linear stability 
code 

Direct numerical 
simulation 

75500 

0 

9.02 

0.4170 

0.4166 

75 500 

74 600 

0 

(1.88, -4048.2) 

(1.96, -4048.2) 

80 700 

0 

1.50 

1.5918 

1.5916 

80700 

73160 

0 

(2.12, -3885.2) 

(2.19, -3885.2) 

80698 

73160 

1.59 

(-26.49, -4844.9) 

(-26.37, -4845.0) 


Additional tests also indicate that the program provides accurate results in the nonlinear regime. 
Another curved channel simulation code written by Zang uses an unsplit procedure (ref. 7) for solving 
the Navier-Stokes equations. When started with identical initial conditions, using a time step of 
A* — 10 -6 , both the split and the unsplit codes indicated that instantaneous complex growth rates 
agreed to about 0.1% at four check points during 400 step simulations. These cases were in the 
nonlinear regime (the maximum streamwise perturbation was 7.2% of the mean centerline velocity), 
hence the instantaneous growth rates were not constant. One further check on the accuracy of 
the split code was obtained for r/ = 0.975, Re = 2466, m = 15, and k — 5. Here a nonlinear, 
nonaxisymmetric wave became periodic with a frequency of 3.095. Using a completely independent 
simulation code, Finlay et al. (ref. 12) obtained a frequency of 3.094. These tests give us confidence 
that our program accurately simulates curved channel flow. 
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4. Comparison of Results 

Here we compare the predictions of the weakly nonlinear theory of DHZ with the results of DNS. 
The number of modes used to resolve the velocity in each direction was chosen to ensure well-resolved 
results. 

In all the simulations the time step (for a full Runge-Kutta step) was At = 10“ 5 . The resultant 
CFL numbers in the simulations were all well below the Runge-Kutta stability limit; most were less 
than half the limit. For cases in which a TS wave was included, approximately 155 full Runge-Kutta 
steps were taken for each TS wave period. 

Sufficient spatial resolution in the simulation can be ensured by using the guideline suggested by 
Krist and Zang (ref. 13). They suggest that “grid refinement is needed in any direction when the tail 
of the energy spectrum reaches 10~ 8 of the low-frequency value.” This guarantees that truncation 
errors in the velocity will be less than 0.01%. Such detailed resolution is not necessary for the 
purposes of these simulations, though unless otherwise specified, the Krist and Zang guideline was 
followed. In what follows, M Chebyshev modes imply that there are M + 1 Chebyshev collocation 
points across the channel; N Fourier modes imply that there are N collocation points in the given 
direction so that the Fourier wave numbers n have a range of —N/ 2 + 1 < nL/( 2n) < N/2 — 1, where 
L is the length of the domain. Typically, 48 Chebyshev modes are used to resolve the flow in the 
wall-normal (radial) direction. In the streamwise (azimuthal) and spanwise (axial) directions, either 
2, 6, or 12 Fourier modes were used, depending on the purpose of the given simulation. The code 
is structured such that at least two modes must be used in each direction. In cases where no TS 
mode is included, two modes are used in the streamwise direction. Similarly, where no Dean mode 
is included, two modes are used in the spanwise direction. 

All cases considered perturbations from the neutrally stable point: Rq = 75000, A = 

A c = 2.179 x 10” 5 , m = 74 600, and k = 9.02. The corresponding dimensional streamwise wave 
number was d ♦ ajs = 1*62. When the curvature ratio was perturbed, the simulations were per- 
formed keeping the value of d • c*ts constant. Qualitatively similar results were obtained at other 
base Reynolds numbers. The cases presented here are considered representative of other parameter 
sets. For the work reported here, the values of the Landau coefficients corresponding to equations (9) 
and (10) are 


01 = 0.908 x 10~ 3 
Si = -0.623 x 10 6 
m = 0.118 x 10 2 
7 ! = 0.155 x 10 7 


fo = 0-355 x 10~ 2 
6 2 = -0.744 x 10 7 
m = -0.851 x 10~ 3 

72 = 0 


Both positive and negative perturbations to the critical curvature ratio were studied. Each is 
presented separately. 

4.1. Positive Curvature Perturbation 

Here the weakly nonlinear theory predicts that TS and Dean modes can attain individual 
equilibrium states, but should they both be present initially, the TS mode will decay and the Dean 
mode will reach its steady state. Unless otherwise specified, the simulations here had e 2 A = 10~ 7 
and Re = Rq -f e 2 R\ = 75 500. In order to keep the streamwise wave number d ♦ c*xs the same as 
that at the critical curvature ratio, equation (8) was used to obtain m = 74 257. This value of m 
was used for all simulations at this curvature ratio. Three separate situations were considered: (1) 
only a Dean disturbance, (2) only a TS wave, and (3) both Dean and TS modes. Note that both 
TS and Dean modes are linearly unstable in this range. 

When only a Dean disturbance is included as part of the initial conditions, the weakly nonlinear 
theory suggests that an equilibrium state ought to be obtained. From equation (9) with the 
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left-hand side equal to zero and \Y\ 2 = 0, we can easily derive the equilibrium state, e 2 \X \ 2 = 
— (e 2 i?l/?l + c 2 A7i)/5i. The corresponding rms value of the streamwise velocity component on the 
centerline is v Yms = e\X\y/2 = 1.41 x 10” 3 . We initiated the simulation with the Dean eigenfunction 
obtained from the linear stability code with an initial strength equal to that predicted by the weakly 
nonlinear theory. Its rms velocity components are given in figure 3. Twelve Fourier modes were 
used in the spanwise direction, though half the modes had energies which were more than 8 orders 
of magnitude less than the primary mode. These higher order modes had a negligible effect on the 
simulation. At the start of the simulation, the instantaneous growth rate of the Dean disturbance 
was approximately equal to its linear value, 0.57. As the simulation progressed, the growth rate 
continually decreased. At time t = 0.13825, the instantaneous growth rate was 0.07 and the 
simulation was stopped with the centerline streamwise rms velocity u' ms — 1.36 x 10 -3 . In this 
case, the difference between the equilibrium state predicted by the weakly nonlinear theory and that 
obtained from the DNS is less than 5% and is considered quite good. 

Initializing this flow with only a TS wave leads to the TS equilibrium state. From equation (10) 
with \X\ = 0 and e 2 |y| 2 = -(e 2 i?i/?2 + e 2 A7 2 )/^2. one finds the rms value of the radial velocity 
component on the centerline is u' ms = e\Y\y/2 = 8.56 x 10“ 4 . The simulation was initiated with 
u'ms = 8.27 x 10~ 4 and was stopped at time t = 0.7465 when the growth rate had slowed to 0.0018. 
At this time, the centerline vertical rms velocity is u' Yms = 9.52 x 10 -4 . Figure 4 illustrates the time 
evolution of the amplitude of the wave. Note that as the amplitude approaches its equilibrium value, 
the growth rate decreases. Over 480 periods of the TS wave were simulated for the data presented in 
figure 4. The equilibrium amplitude predicted by the weakly nonlinear theory is shown by the dashed 
line. The difference between the equilibrium amplitude predicted by the weakly nonlinear theory and 
that obtained by the DNS is about 15%, which is considered acceptable agreement by Philip Hall. 
The 12 Fourier modes used in the streamwise direction were sufficient for an 8-order-of-magnitude 
drop in the energy. 

We performed several simulations with both Dean and TS perturbations included in the initial 
conditions. In these, we typically used 6 modes in both the streamwise and the spanwise directions. 
As usual, 48 Chebyshev modes were used in the wall-normal direction. In these simulations, the 
spanwise resolution was only sufficient to allow an energy drop of 10 5 , while the streamwise resolution 
allowed as little as a 10 4 drop in the energy. This would lead us to suspect errors in the velocity on 
the order of 1%. To determine whether this could lead to erroneous physical results, we resimulated 
our case which had the largest velocity perturbations, using 12 modes in the streamwise and spanwise 
directions and 64 modes in the wall-normal direction. The energy spectrum in the highly resolved 
case exhibited an energy drop of 10 9 in the streamwise direction and 10 14 in the spanwise direction. 
The evolution of the disturbance amplitudes was the same as that observed in the less resolved case. 

In the remaining parts of this subsection, unless otherwise specified, we initiated the Dean 
disturbances with v Yms = 1.41 x 10“ 3 on the centerline. This value is slightly more than its 
equilibrium amplitude which we obtained above. The weakly nonlinear theory predicts that for 
the given conditions, the Dean disturbance will evolve towards this equilibrium state. 

In the first case, the initial strength of the TS disturbance was = 6.93 x 10 -5 , approximately 
1/12 its predicted equilibrium value. The long time behavior of this case required 12 Fourier modes 
in the spanwise direction. The amplitude evolution of the disturbances is illustrated in figure 5(a). 
Here we plot the amplitude of both the Dean and TS disturbances, normalized by their respective 
amplitudes at time t = 0. The Dean disturbance develops towards its equilibrium state while the 
TS disturbance decays. This behavior is exactly what one would expect from the weakly nonlinear 
theory. 

In the second case, the initial strength of the TS disturbance was u' Ttns = 6.93 x 10 -4 , 
approximately 5/6 its predicted equilibrium value. Figure 5(b) shows the amplitude histories for 
this case. Unlike the previous situation, here the TS wave grows slightly and the Dean disturbance 
experiences rapid decay. The decay of the Dean disturbance is not predicted by the weakly nonlinear 
theory. In an attempt to understand this discrepancy we look at the rms distribution of velocities 
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in the Fourier mode which characterizes the Dean disturbance in figure 6. Note that the vertical 
component of velocity has two local maxima in its distribution, and the spanwise component has 
three local maxima. Figure 3 illustrates the velocity components at time t = 0. Here the vertical 
and spanwise velocities have one and two local maxima, respectively. At the end of the simulation, 
the TS rms velocity distributions were the same as they were initially. In the previous case with the 
weaker TS wave, the Dean disturbance maintained its original structure throughout the simulation. 
The stronger TS disturbance in the present simulation distorts the cross-stream velocities in the 
Dean mode. Without the appropriate cross-stream velocity distribution, the streamwise velocity 
component of the Dean mode decays. Note that this happens quickly; the data in figure 5(b) 
represent approximately 15 TS wave periods. We continued this simulation with 12 modes in the 
streamwise direction. As seen in figure 7, the flow develops into a state with the TS wave at its 
equilibrium amplitude. Another case including both TS and Dean modes of comparable amplitude, 
but with e 2 i?i = 100 and e 2 A — 10 -8 , evolved similarly. 

Another simulation was performed with precisely the same initial conditions as were used in 
the above paragraph. Instead of allowing the Dean disturbance to evolve naturally in time, we 
artificially forced its original amplitude and velocity distribution to be constant. Figure 8 shows 
the development of both the Dean and the TS disturbances. Unlike the previous case, here the TS 
disturbance decays in time, as one would expect from the weakly nonlinear theory. This indicates 
that the decay of the TS mode is directly linked to the presence of a strong, self-sustaining Dean 
vortex. 

One more case of interest is when the Reynolds number perturbation is 2000, i.e., Re = 
Ro 4- e 2 R\ = 77000. The perturbation from the critical curvature ratio remained e 2 A = 10“ 7 . 
Initially the amplitudes of the Dean and TS disturbances are such that for their respective 
disturbances, v f rms = 1.38 x 10 -3 and u' rms = 9.70 x 10” 4 on the centerline. Figure 9 illustrates the 
energy in the Dean mode, the TS mode, and the oblique mode with wave vector (m = 74257, k = 
9.02). Here energy is plotted instead of the magnitude of velocity at a single point in the channel. 
Since this case initially evolves quite similarly to the case above in which strong Dean and TS 
disturbances were allowed to develop, the velocity distribution in the Dean mode changes its structure 
in time, so the integrated energy in each mode is a much better measure of its size. Before t zz 0.08, 
the drastic changes in the slope of the energy-versus-time plot for both the Dean and the oblique 
modes are correlated with large changes in the relative magnitudes of the streamwise and wall- 
normal components of the velocity in the Dean mode. The reason(s) for these changes is (are) 
not understood. For t > 0.08, both the oblique and the Dean modes become more energetic at 
approximately the same rate. This is due to a secondary instability similar to that described by 
Herbert (refs. 14 and 15) for the plane channel. According to the secondary instability theory, the 
oblique mode and the mode corresponding to the Dean vortex become unstable when the amplitude of 
TS disturbance exceeds some critical value that depends on the Reynolds number, the spanwise wave 
number, and the geometry of the flow. As can be seen from our results, this secondary instability 
is quite strong and will likely lead to the development of turbulence in the flow. Unlike the case 
above in which the TS mode achieved an equilibrium state before its amplitude was great enough 
to support a secondary instability for the spanwise wave numbers simulated, no simple equilibrium 
or periodic states are likely to develop in this flow. 

An additional simulation was performed to investigate the possibility that higher order modes 
which are included in the simulation but not in the weakly nonlinear theory could be responsible for 
the discrepancy. The theory includes contributions from the Dean fundamental and first harmonic, 
the TS fundamental and first harmonic, the mean flow distortion, and the interaction between the 
Dean fundamental and the TS fundamental. These modes were allowed to evolve as usual. The 
modes corresponding to the interaction of either harmonic with either fundamental or each other 
were artificially set to zero at every time step. In this way, only the modes which the theory 
explicitly considers were permitted to evolve in the simulation. We initialized a simulation with 
the TS disturbance approximately equal to its equilibrium state and obtained results which were 
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qualitatively similar to what we had without suppressing any modes; the Dean mode decayed rapidly 
and the TS mode remained approximately constant. 

4.2. Negative Curvature Perturbation 

We explored one case in which the perturbation to the channel curvature was the opposite of that 
used above. With e 2 A = —1.0 x 10 -7 and Re = 75 100, the weakly nonlinear theory predicts that the 
combined TS-Dean mode state dominates. As before, the value of m needed to be adjusted to keep 
d- ax S constant. Using equation (8), we found m = 74 942. Using these parameters, the linear growth 
rate of the TS wave is 0.33, while that of the Dean mode is —0.10. Weakly nonlinear theory predicts 
a combined equilibrium state with the Dean t>£ ms = 3.08 x 10 -4 and the TS u' ms = 1.68 x 10 -6 . 

The simulation was done using six modes in both the streamwise and the spanwise directions. 
The energy spectrum showed a drop of 10 12 in the streamwise direction and 10 6 in the spanwise 
direction. In view of our experience with the positive curvature perturbation cases, this resolution 
was deemed sufficient for the purpose of determining the existence of a combined mode state. The 
amplitudes of the respective modes, divided by their initial amplitudes, are illustrated in figure 10. 
The Dean disturbance decays at the rate of -0.11, slightly faster than its linear rate obtained from 
the finite-difference code. The growth rate of the TS disturbance is 0.22, slightly less than its linear 
rate. The DNS does not indicate that a combined mode state exists with these parameters. It is 
possible that such a state exists with amplitudes of the respective disturbances somewhat different 
from those predicted by DHZ; however, the computational cost of finding such a state without 
additional information would be excessively high. 

5. Discussion 

Comparison of the results of the weakly nonlinear theory of DHZ with the results of DNS indicates 
that the theory and the simulation disagree with respect to some aspects of the flow. 

When the curvature parameter perturbation of the channel is positive, any small, positive 
Reynolds number perturbation allows equilibrium states for isolated Dean and TS modes. The 
amplitudes of these states obtained with the weakly nonlinear theory and the DNS agree to 
acceptable accuracy. When both modes are included as part of the initial conditions, the weakly 
nonlinear theory predicts that the ultimate state will contain only a Dean mode in equilibrium; 
the TS mode is predicted to decay to zero. The DNS results indicate that this will occur if the 
initial Dean amplitude approximates its equilibrium amplitude and the TS amplitude is small. 
When the TS amplitude is close to its equilibrium state, the cross-stream velocity perturbations 
in the Dean mode lose their coherence and the Dean amplitude decays. If the structure and 
amplitude of the Dean mode are artificially maintained in their original condition, the TS mode 
will decay even though its initial amplitude is not small. We suggest two possible sources of 
this discrepancy: (1) the cases simulated were out of the range of validity of the theory, or 

(2) one of the Landau coefficients in equations (9) and (10) is incorrect. We discuss each of these 
possibilities below. 

All weakly nonlinear theories suffer uncertainty as to their range of validity. There is no a priori 
way of knowing this range; it can only be checked a posteriori by comparison with some higher order 
theoretical predictions, physical experiments, or computations. The DHZ theory perturbs both the 
Reynolds number and the curvature ratio. In the cases presented above, the Reynolds number 
perturbation was 500 from a base Reynolds number of 75 000; the curvature parameter perturbation 
was 0.01 x 10 -5 from a critical value of 2.179 x 10 -5 . The Dean mode also decayed for a case 
with much smaller Reynolds number and curvature ratio perturbations. While it is possible for the 
perturbations used to be outside the valid range of the theory, such a limited range of validity would 
give the theory little predictive or explanatory value. 

A numerical error in the calculation of one of the Landau coefficients in DHZ could explain the 
discrepancy between the theory and the DNS. Since equilibrium amplitudes of both Dean and TS 
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disturbances agree with the two approaches, it seems likely that the coefficients ft, 7i, 61, 72 > and 

r/2 are approximately correct. When the Dean mode maintains itself (either artificially or because 
the TS mode is too weak to have a large impact on it), both the theory and the DNS agree that 
TS disturbances decay; hence, the coefficient responsible for the effect of the Dean mode on the TS 
mode, 62, seems to be correct. The DNS and the weakly nonlinear theory disagree on the effect 
of the TS mode on the Dean mode. This is felt only through 7/1, hence it is possible that 7/1 was 
computed incorrectly. 

A numerical error in the computation of 7/1 would also help explain the discrepancy between the 
predictions of the weakly nonlinear theory and the DNS for the case of negative curvature ratio 
perturbations. In the parameter range simulated, the weakly nonlinear theory predicts that both 
TS and Dean modes should coexist. The Dean mode by itself decays; there is no Dean equilibrium 
state. Only the influence of a TS wave, felt through the coefficient r )\ , could make the sustained 
existence of a Dean mode possible. The DNS results do not indicate that the two modes coexist 
in equilibrium for the set of parameters studied. As before, the value of 771 seems to be behind the 
discrepancy. 

6. Conclusions 

We performed direct numerical simulations (DNS) of curved channel flow in order to validate the 
weakly nonlinear theory of Daudpota, Hall, and Zang (DHZ). The two approaches provide amplitudes 
of equilibrium states for individual Dean and Tollmien-Schlichting (TS) disturbances which agree 
to acceptable accuracy. Unfortunately, the theory and the DNS are in conflict as to the effects of 
interaction between the two modes. Two possible sources of the disagreement have been discussed. 
It is possible that the weakly nonlinear theory is valid only for a very narrow range of parameters and 
that its use here is out of that range. A more likely possibility is indicated by the character of the 
discrepancies, which suggests that the effect of the TS disturbance on the Dean mode is not being 
correctly predicted by the theory. The reverse interaction, that of the Dean on the TS, does seem 
to be correctly predicted; hence, we believe that while the bulk of the theory was done correctly, the 
Landau coefficient, which reflects the influence of the TS on the Dean, was incorrectly calculated by 
DHZ. 


NASA Langley Research Center 
Hampton, VA 23665-5225 
June 27, 1989 
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(a) TS disturbances. 



(b) Dean disturbances. 
Figure 2. Neutral stability curves. 
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(a) Initial TS amplitude is approximately 1/12 its predicted equilibrium amplitude. 



(b) Initial TS amplitude is approximately 5/6 its predicted equilibrium amplitude. 

Figure 5. Time history of amplitude of Dean and TS modes, where A is the amplitude normalized by the 
respective initial amplitude. 
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Figure 6. Rms velocity distribution in the Dean mode. Re = 75 500; e 2 A = 10 7 ; t = 0.021; initial TS amplitude 
is approximately 5/6 its predicted equilibrium amplitude. 
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Figure 8. Time history of amplitude of Dean and TS modes, where A is the amplitude normalized by the 
respective initial amplitudes. Initial TS amplitude is approximately 5/6 its predicted equilibrium amplitude; 
the Dean disturbance is held constant in time. 
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Figure 9. Time history of energy in Dean, TS, and oblique modes. Re — 77 000. 



Figure 10. Time history of amplitude of Dean and TS modes, where A is the amplitude normalized by the 
respective initial amplitudes. Re = 75 100; e 2 A = —1.0 x 10 -7 ; weakly nonlinear theory predicts a combined 
mode state. 
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